Storm Track Analysis¶
In [32]:
#Create tracking info ensemble
UM133_ens = xr.concat(UM133.dataset.values(), dim='ens')
UM044_ens = xr.concat(UM044.dataset.values(), dim='ens')
CPOL = OBS.dataset
In [33]:
#Define some functions
def get_quant(df, **kwargs):
try:
quant = kwargs['quant']
except KeyError:
quant={1:(0,0.2), 2:(0.2,0.4), 3:(0.4, 0.6), 4:(0.6,0.8), 5:(0.8,1)}
try:
rain = kwargs['rain']
except KeyError:
rain='avg_mean'
for n, q in quant.items():
Q = df[rain].quantile(q)
df.quant.loc[(df[rain]>= Q[q[0]]) & (df[rain]<Q[q[1]]) ] = int(n)
return df
def ravel(ens_data, run='Cpol', **kwargs):
pdf = {}
for var in dict(ens_data.variables).keys():
if not var.startswith('dim') :
pdf[var] = ens_data[var].values.ravel()
pdf[var] = pdf[var][pdf[var] != np.nan]
df = pd.DataFrame(pdf).dropna()
run_n = [run for i in range(len(df))]
df['run'] = run_n
df['quant'] = np.zeros(len(df), dtype='i')
df.index = np.arange(len(df))+1
return get_quant(df, **kwargs)
UM133_pdf, UM044_pdf, CPOL_pdf = ravel(UM133_ens, 'UM 1.33km'), ravel(UM044_ens, 'UM 0.44km'), ravel(CPOL, 'Cpol')
def matplotlib_to_plotly(cmap, pl_entries, rgb=True):
h = 1.0/(pl_entries-1)
pl_colorscale = []
for k in range(pl_entries):
C = list(map(np.uint8, np.array(cmap(k*h)[:3])*255))
if rgb:
pl_colorscale.append([k*h, 'rgb'+str((C[0], C[1], C[2]))])
else:
pl_colorscale.append(tuple(np.array(C)/255))
return pl_colorscale
The tracking algorithm is a fork of the Tint (Tint is not Titan) tracking algorithm (http://openradarscience.org/TINT/). The framework has been modified that it can be also applied to model output data that is not not stored in Py-ART radar data container.
- The analysis is based averages of one week worth of strome cell tracking:
1. Median of strom area, duration, avg., max. rain-rates and # of storm cells¶
In [34]:
#Plot the Medians
um044_m =np.array(list(UM044_t.num.values())).mean(dtype='i')
um133_m =np.array(list(UM133_t.num.values())).mean(dtype='i')
cpol_m =np.array(list(CPOL_t.num.values())).mean(dtype='i')
var=['avg_area','dur', 'avg_mean', 'max_mean']
names=['area', 'duration', 'avg-rain', 'max-rain', '# cells']
#print(CPOL_pdf[var].median())
medians = pd.DataFrame({'a': list(CPOL_pdf[var].median().values)+[int(cpol_m)],
'b': list(UM044_pdf[var].median().values)+[int(um044_m)],
'c': list(UM133_pdf[var].median().values)+[int(um133_m)]} )
#index=('Area','Duration', 'Mean-Rain', 'Max-Rain'))
medians.index=names
medians.columns=['Cpol', 'UM 1.33km', 'UM 0.44km']
#print('Medians:')
medians.round(2)
Out[34]:
2. 2D-Histograms (storm area vs duration):¶
In [46]:
#Create Hex-Bin plot
mpld3.disable_notebook()
var=['avg_area','dur', 'max_mean', 'avg_mean']
medians = pd.DataFrame({'CPOL':CPOL_pdf[var].median(),
'UM 1.33km':UM133_pdf[var].median(),
'UM 0.44km':UM044_pdf[var].median()})
#medians.loc['avg_area'] /= 2.5**2
histdata = [UM044_pdf[var].dropna(), UM133_pdf[var].dropna(), CPOL_pdf[var].dropna()][::-1]
titles = ['UM 0.44km ens', 'UM 1.33km ens', 'CPOL'][::-1]
fig = plt.figure(figsize=(10,8))
colm = colmap2.Blues
colm.set_under('w', alpha=0)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True)
fig.subplots_adjust(bottom=0.07, right=0.9, top=0.35, wspace=0.01)
cbar_ax = fig.add_axes([0.15, 0.02, 0.7, 0.01])
hexbin_data = []
gridsize = 10
vmin=0.0001
vmax=0.1
nticks=10
YMax, XMax = 260, 70*2.5**2
for i, ax in enumerate((ax1, ax2, ax3)):
data = histdata[i][var[:2]]
#data[var[0]] /= 2.5**2
#data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
X = data[var[0]].values
Y = data[var[1]].values
ax.set_ylim(0,YMax)
ax.set_xlim(0,XMax)
hb = ax.hexbin(X, Y, gridsize=gridsize, extent=(0,XMax,0,YMax))
hexbin_data.append(np.ones_like(Y, dtype=np.float) / hb.get_array().sum())
plt.cla()
medians = OrderedDict()
for i, ax in enumerate((ax1, ax2, ax3)):
data = histdata[i][var[:2]]
#data[var[0]] /= 2.5**2
#data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
X = data[var[0]].values
Y = data[var[1]].values
ax.set_ylim(0,YMax)
ax.set_xlim(0,XMax)
im = ax.hexbin(X, Y, gridsize=gridsize, cmap=colm, marginals=False, extent=(0,XMax,0,YMax),
vmin=vmin, vmax=vmax, C=hexbin_data[i], reduce_C_function=np.sum)
ax.set_title(titles[i], fontsize=24)
ax.grid(color='w', linestyle='', linewidth=0)
ax.tick_params(labelsize=24)
ax.xaxis.set_ticks(ax.xaxis.get_ticklocs()[:-1])
x, y = histdata[i][var[0]].median(), histdata[i][var[1]].median()
z, zz= histdata[i][var[2]].median(), histdata[i][var[-1]].median()
sx, sy = histdata[i][var[0]].std(), histdata[i][var[1]].std()
medians[titles[i]] = '%2.1f km^2, %2i min (max: %2.1f mm/h, mean: %2.1f mm/h);'%(x,y, z, zz)
ax.hlines(y*np.ones_like(histdata[i][var[0]]),0,x, 'firebrick', lw=3)
ax.vlines(x*np.ones_like(histdata[i][var[1]]),0,y, 'firebrick', lw=3)
ax.scatter([x], [y], marker='o', s=400, c='firebrick', alpha=0.8)
if i == 0:
ax.set_ylabel('Duration [min]', fontsize=24)
ax.set_xlabel('Rainfall Area [km$^2$]', fontsize=24)
ary=im.get_array()/im.get_array().sum() * 100
im.set_array = ary
cbar = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Density [ ]',size=24)
#cbar.set_ticks(np.linspace(vmin, vmax, nticks).round(2))
#cbar.set_ticklabels(np.linspace(vmin, vmax, nticks).round(2))
#print('Medians:')
#medians
#print(' '.join(['%s: %s'%(k, v) for (k,v) in medians.items()]))
plt.show()
In [47]:
colm = matplotlib_to_plotly(CubeYF_20.get_mpl_colormap(N=8, gamma=2.0),8, rgb=False)
3. Storm tracks above the 80th percentiles¶
- marks start of track, colors for different ensemble member
In [42]:
#Plot the tracks
mpld3.enable_notebook()
fig = plt.figure(figsize=(13,8))
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.9,left=0.01, hspace=0, wspace=0)
#cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
tp = 'mean'
num=80
titels = ['CPOL', 'UM 1.33km', 'UM 0.44km']
m = None
for i, tracks in enumerate((CPOL_t, UM133_t, UM044_t)):
ax = fig.add_subplot(1,3,i+1)
ax.set_title(titels[i], fontsize=18)
for ii, tr in enumerate(tracks.dataset.keys()):
if ii == 0:
draw_map = None
else:
draw_map = m
perc = tracks.percentiles[tr][tp][num]
ax, m, im = plot_traj(tracks.dataset[tr], lon, lat, ax=ax, mintrace=2, thresh_val=perc, thresh='mean',
color=colm[ii], draw_map=draw_map, basemap_res='f', lw=0.5, size=20)
plt.show()
4. Classification of strom area, rain-rate and duration by rain-rate quintiles¶
In [60]:
mpld3.disable_notebook()
variables=dict(dur=('Duration [min]',300), avg_area=('Area [km$^2$]', 100*2.5**2), avg_mean=('Rain-rate [mm/h]',12))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
label, ylim = prop
npl = len(list(variables.keys()))
i += 1
ax = fig.add_subplot(npl,1,i)
ax = sns.boxplot(x="quant", y=y, hue="run", data=df_stack, palette="muted", ax=ax)
#ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
ax.set_xlim(0.5,5.5)
ax.set_ylim(0,ylim)
ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
ax.set_xlabel('Quintiles', fontsize=26)
ax.set_ylabel(label, fontsize=26)
5. Stormtrack Percentiles vs Rain-rate¶
In [58]:
P = list(np.arange(0,99))+[99, 99.9, 99.99, 99.999, 100]
PERC = pd.DataFrame({'Obs':np.percentile(CPOL_pdf['avg_mean'].values, P),
'UM 1.33km': np.percentile(UM133_pdf['avg_mean'].values, P),
'UM 0.44km': np.percentile(UM044_pdf['avg_mean'].values, P)},index=P)
from mpl_toolkits.axes_grid.inset_locator import inset_axes
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
ax.grid(color='r', linestyle='-', linewidth=0)
